function f=find_S_w(w,P_i,W_i,h_w,BB,theta_hat,k_type,tau,packobj)

N=size(W_i,1);

f=NaN*w;
S_w=f;

	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	% UNPACK <packobj>
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
sigma1=packobj.sigma1;
sigma2=packobj.sigma2;
fdist1=packobj.fdist1;
fdist2=packobj.fdist2;
sigma_berkson=packobj.sigma_berkson;
alpha1=packobj.alpha1;
alpha2=packobj.alpha2;
q=packobj.Q_i;
K1=packobj.K1;
K2=packobj.K2;
p_lb=packobj.p_lb;
p_ub=packobj.p_ub;
q_lb=packobj.q_lb;
q_ub=packobj.q_ub;
mu1=0;
mu2=NaN;

	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	% COMPUTE G-inverse
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
		
Ginv_hat=BB*theta_hat;
Ginv_hat_term=(Ginv_hat-tau<=0);

Ginv_hat_term_withB=NaN*Ginv_hat_term;
for ii=1:N

		f_int=@(eps) ((construct_chebyshev_berkson_Ponly_twodim_selBB(P_i(ii)+eps,q(ii), ...
				K1,K2,p_lb,q_lb,p_ub,q_ub,NaN,NaN,mu1,sigma1,fdist1,alpha1)'*theta_hat)<=tau)' .* normpdf(eps,0,sigma_berkson);

		cutoff=4*sigma_berkson; 

		lb=max(-P_i(ii),-cutoff); 
		ub=min(-P_i(ii)+1,+cutoff);
		Ginv_hat_term_withB(ii,:)=integral(f_int,lb,ub,'AbsTol',0,'RelTol',1e-3);
end
	
for jj=1:size(w(:),1)
	K_w=kfunc((W_i*ones(1,length(w(jj)))-ones(N,1)*w(jj))/h_w,k_type);
	S_w(jj) = (1/(h_w))*mean((Ginv_hat_term_withB).*K_w,1);
end

f=S_w;

	%
